###############################
# Analysis of conjoint data 
###############################

rm(list=ls())

# load packages 
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library(Hmisc)

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##############################
# Prepare data
##############################

# load data
d = read.dta("~/Dropbox/Cerrillos 2017/08_replication/conjoint_cerrillos_LAPS.dta")
table(d$failcon)
names(d)
dim(d)

# check data
d$idnum = as.numeric(d$idnum)
table(d$enumerator)
table(d$selected, exclude = NULL)

# Vote past election
table(d$a11)
d$votepastelection = 0
d$votepastelection[d$a11==1]=1
table(d$votepastelection)
round(prop.table(table(d$votepastelection)),2)

# Vote next election
d$a10 = as.numeric(d$a10)
table(d$a10, exclude = NULL)
d$votenextelection = 0
d$votenextelection[d$a10<7]=1
table(d$votenextelection)
round(prop.table(table(d$votenextelection)),2)

d$novotenextelection = 1 - d$votenextelection
table(d$novotenextelection)
round(prop.table(table(d$novotenextelection)),2)

# Name last election
table(d$a12)
d$namepastelection = 0
d$namepastelection[d$a12==1]=1
d$namepastelection[d$a12==2]=1
d$namepastelection[d$a12==3]=1
d$namepastelection[d$a12==4]=1
d$namepastelection[d$a12==5]=1
table(d$namepastelection)
prop.table(table(d$namepastelection))

# Interest in politics
table(d$a3, exclude = NULL)
d$interest = 0
d$interest[d$a3<4]=1
table(d$interest)
prop.table(table(d$interest))

# Likely voter 1
d$likely1 = 0
d$likely1[d$interest==1 & d$votepastelection==1 & d$votenextelection==1]=1
prop.table(table(d$likely1))
d$nolikely1 = 1 - d$likely1

# drop failcon and non-responses
d = subset(d, d$failcon!=1)
d = subset(d, d$selected!=99)
d = subset(d, d$selected!=88)
table(d$selected, exclude = NULL)

# prepare centrist and nonresponse sample
d = d[d$a5>4 & d$a5<7|d$a5==88|d$a5==99,]

# subgroups
d$rightwingcentrist = 0
d$rightwingcentrist[d$a4==2 & d$d2==1 & d$c2==2]=1
table(d$rightwingcentrist, exclude = NULL)

d$leftwingcentrist = 0
d$leftwingcentrist[d$a4==1]=1
d$leftwingcentrist[d$d2==2 & d$c2==2]=1
table(d$leftwingcentrist, exclude = NULL)

d$nocandidates = 0
d$nocandidates[d$rightwingcentrist==0 & d$leftwingcentrist==0]=1
table(d$nocandidates, exclude = NULL)

# subsets 
d_right = d[d$rightwingcentrist==1,]
d_left = d[d$leftwingcentrist==1,]
d_nocandidate = d[d$nocandidates==1,]

# ideological centrists and non-responses
as.numeric((dim(d_right)[1]) + (dim(d_left)[1]))/as.numeric(dim(d)[1])

# check likely
table(d_left$likely1)
table(d_right$likely1)
table(d_nocandidate$likely1)

###################
# Main results
###################

# Left-wing candidate
describe(d_left$idnum)
left <- lm(d_left$outcome ~ d_left$atideology + d_left$atprofession + d_left$atage)
left_c <- round(coeftest(left, vcov = vcovCluster(left, cluster = d_left$idnum)),4)
left_c

# Right-wing candidate
describe(d_right$idnum)
right <- lm(d_right$outcome ~ d_right$atideology + d_right$atprofession + d_right$atage)
right_c <- round(coeftest(right, vcov = vcovCluster(right, cluster = d_right$idnum)),4)
right_c

# No candidate
describe(d_nocandidate$idnum)
nocandidate <- lm(d_nocandidate$outcome ~ d_nocandidate$atideology + d_nocandidate$atprofession + d_nocandidate$atage)
nocandidate_c <- round(coeftest(nocandidate, vcov = vcovCluster(nocandidate, cluster = d_nocandidate$idnum)),4)
nocandidate_c

#####################
# Prepare plots
#####################

# left
left_c # regression results
coeff = left_c [1:6,1] # save coefficientes
coeff2 = cbind(coeff) # matrix of coefficientes
c_se = left_c [1:6,2] # save cluster standard errors
c_se2 = cbind(c_se) # matrix of cluster standard errors
results = data.frame(coeff2,c_se2) # generate dataset 
colnames(results) = c("y1","r1")  # change column names
results = results[-1,] # drop intercept 
baseline =  c(0,0) # gen baseline rows
results <- rbind(baseline,results[1, ], baseline,results[2:3, ],baseline,results[4:5, ]) # add baselines
rownames(results) = c("1b.atideology","2.atideology",
                      "1b.atprofession","2.atprofession","3.atprofession",
                      "1b.atage","2.atage","3.atage")
results 
write.table(results, "~/Dropbox/Cerrillos 2017/08_replication/013_centrist_left.txt", sep="\t")

# right
right_c # regression results
coeff = right_c [1:6,1] # save coefficientes
coeff2 = cbind(coeff) # matrix of coefficientes
c_se = right_c [1:6,2] # save cluster standard errors
c_se2 = cbind(c_se) # matrix of cluster standard errors
results = data.frame(coeff2,c_se2) # generate dataset 
colnames(results) = c("y1","r1")  # change column names
results = results[-1,] # drop intercept 
baseline =  c(0,0) # gen baseline rows
results <- rbind(baseline,results[1, ], baseline,results[2:3, ],baseline,results[4:5, ]) # add baselines
rownames(results) = c("1b.atideology","2.atideology",
                      "1b.atprofession","2.atprofession","3.atprofession",
                      "1b.atage","2.atage","3.atage")
results 
write.table(results, "~/Dropbox/Cerrillos 2017/08_replication/014_centrist_right.txt", sep="\t")

# nocandidate
nocandidate_c # regression results
coeff = nocandidate_c [1:6,1] # save coefficientes
coeff2 = cbind(coeff) # matrix of coefficientes
c_se = nocandidate_c [1:6,2] # save cluster standard errors
c_se2 = cbind(c_se) # matrix of cluster standard errors
results = data.frame(coeff2,c_se2) # generate dataset 
colnames(results) = c("y1","r1")  # change column names
results = results[-1,] # drop intercept 
baseline =  c(0,0) # gen baseline rows
results <- rbind(baseline,results[1, ], baseline,results[2:3, ],baseline,results[4:5, ]) # add baselines
rownames(results) = c("1b.atideology","2.atideology",
                      "1b.atprofession","2.atprofession","3.atprofession",
                      "1b.atage","2.atage","3.atage")
results 
write.table(results, "~/Dropbox/Cerrillos 2017/08_replication/015_centrist_nocandidate.txt", sep="\t")

